/* Modular exponentiation using repeated squaring . */
/* That is, we want to compute a^b mod n. */
modexp(a, b, n) = { \
    local(d, bin); \
    d = 1; \
    bin = binary(b); \
    for (i = 1, length(bin), \
        d = Mod(d*d, n); \
        if (bin[i] == 1, \
            d = Mod(d*a, n); \
        ); \
    ); \
    return(d); \
}

modexp2(b, e, m) = {
	local(r, bin);
	r = 1;
	bin = binary (e);
	for ( i = 1, length(bin),
		r = (r*r) % m;
		if( bin[i] == 1,
			r = (r*b) % m;
		);
	);
	return (r);
}
